Non-Intestinal Microbial Signatures within Stool as Predictors of Cancer Immunotherapy Outcome

Introduction

Despite the increasing use of immunotherapy in treating advanced malignancies, clinical outcomes remain heterogeneous, with many patients exhibiting intrinsic resistance or eventual disease recurrence following immune checkpoint inhibitor (ICI) therapy [Robert et al., 2015]. This highlights the critical need for predictive biomarkers to optimize treatment selection and personalize therapeutic strategies. Accumulating evidence indicates a strong correlation between gut microbial composition and responsiveness to cancer immunotherapy [Peng et al., 2020]. While advancements in targeted therapies and immunotherapy have improved survival outcomes [Switzer et al., 2022], ICIs targeting cytotoxic T-lymphocyte-associated antigen 4 (CTLA-4) and programmed cell death protein 1 (PD-1) demonstrate limited efficacy in approximately 50% of patients with metastatic disease [Larkin et al., 2015] and carry the risk of immune-related adverse events [Robert et al., 2019].

Consequently, elucidating the mechanistic relationship between the gut microbiome and antitumor immunity has become a priority in translational cancer research. Indeed, the gut microbiome is increasingly recognized as a critical modulator of anti-tumor immunity and, consequently, response to cancer immunotherapy. Initial investigations, utilizing both preclinical models [Sivan et al., 2015] and clinical cohorts, particularly those with melanoma [Routy et al., 2018], have demonstrated a compelling link between gut microbial composition and immunotherapy efficacy. This link has been further supported by fecal microbiota transplantation studies in both gnotobiotic mice [Gopalakrishnan et al., 2018] and human patients [Baruch et al., 2021]. While much of the early research focused on melanoma [Tsakmaklis et al., 2023], current efforts are shifting towards identifying conserved microbial biomarkers and mechanisms applicable across diverse cancers [Gunjur et al., 2024; Lin et al., 2025], predicated on the hypothesis that fundamental interactions between the gut microbiome and the immune system are broadly relevant to immunotherapy response.

Despite numerous published studies, the biological mechanisms by which gut microbes regulate the immune system in the context of melanoma immunotherapy remain incompletely understood. A recent meta-analysis identified robust stool metagenomic markers associated with successful melanoma immunotherapy [Olekhnovich et al., 2023]. Notably, no taxonomic predictors of negative response were identified. Expanding on our earlier investigations, we employed genome-resolved metagenomics to gain a deeper understanding of the gut microbiota’s role in modulating immunotherapy response in malignant tumors. This analysis revealed a link between systemic metabolic depletion in the gut microbiome and decreased efficacy of immunotherapy in patients with metastatic melanoma [Zakharevich et al., 2024].

We expanded our application of genome-resolved metagenomics, analyzing 624 stool metagenomes from 11 published studies encompassing over five cancer types and representing diverse geographical regions across three continents. This analysis identified robust markers capable of differentiating patients based on immunotherapy response. These markers were validated using statistical analysis and machine learning techniques focused on logarithmic ratios of identified features. Notably, the analysis revealed non-intestinal traits such as oral- and food-derived bacteria within the marker list predominantly associated with negative treatment outcomes, potentially reflecting the previously observed depletion of gut microbial load associated with oral microbiota enrichment [Liao et al., 2024]. These findings underscore the importance of considering the metabolic landscape of the gut microbiome – and the surprising influence of non-intestinal microbial sources – when predicting and optimizing immunotherapy response, paving the way for the development of novel therapeutic strategies targeting microbiome-immune interactions.

References:

Materials and methods

Data acquisition and preprocessing

All 624 stool metagenomic samples from 11 published studies were downloaded from the NCBI/EBI databases using Kingfisher v0.4.1 [Woodcroft, 2022]. Data quality was assessed with FastQC v0.12.1 [Andrews et al., 2023], and raw reads were processed with fastp v0.23.4 to remove low-quality sequences [Chen et al., 2018]. Human-derived sequences were filtered out using HISAT2 v2.2 [Kim et al., 2019] with the GRCh37 human genome (Release 47) as a reference [Schneider et al, 2017].

Assembly and binning

Metagenomic assembly was performed using MEGAHIT v1.2.9 [Li et al., 2015], retaining contigs longer than 1000 bp. Filtered contigs were aligned to metagenomic reads using HISAT2 v2.2 [Kim et al., 2019]. Binning was conducted in two stages: (1) initial binning with MetaBAT2 v2.12.1 [Kang et al., 2019], MaxBin2 v2.2.7 [Wu et al., 2016], and Semibin2 v2.1.0 [Pan et al., 2023], followed by (2) refinement with DAS Tool v1.1.7 to improve bin quality [Sieber et al., 2018]. Bins were quality-checked with CheckM v1.2.3 [Parks et al., 2015] and dereplicated at 98% nucleotide identity using dRep v3.4.5 [Olm et al., 2017] to generate operational genomic units (OGUs).

Taxonomic profiling and analysis

Taxonomic annotation of bins was performed using GTDB-Tk v2.1.1 [Chaumeil et al., 2022] with the GTDB r207 database [Parks et al., 2022]. Additional data processing utilized Samtools v1.17 [Li et al., 2009], BEDTools v2.31.0 [Quinlan et al., 2010], and BBMap v39.06 [Bushnell et al., 2014]. Finally, OGU abundance profiles were generated by mapping reads to the OGU catalog using HISAT2 v2.2, followed processong with InStrain v1.9.0 [Olm et al., 2021]

Marker OGU discovery

The identification of OGUs associated with immunotherapy outcomes followed an established analytical framework from our previous works [Olekhnovich et al., 2023; Zakharevich et al., 2024]. Differential rankings was first performed using Songbird [Morton et al., 2019] to identify OGUs showing relative abundance variations between experimental groups, applying a conservative absolute differential value threshold of > 0.3. For candidate OGUs meeting this criterion, we subsequently calculated log-ratio abundances using Qurro [Fedarko et al., 2020] and determined statistical significance through Wilcoxon rank-sum tests implemented in the R statistical environment. The biomarker selection process incorporated stringent cross-validation criteria to ensure robust identification of microbial signatures. OGUs demonstrating consistent positive associations with therapeutic response across multiple datasets were retained as potential beneficial biomarkers, while any evidence of negative association with treatment outcome in any dataset resulted in automatic exclusion regardless of other positive associations.This approach enabled simultaneous identification of two clinically meaningful biomarker categories: microbial taxa positively correlated with successful immunotherapy outcomes and those associated with adverse therapeutic responses. The methodology emphasizes reproducibility through multi-dataset validation and maintains rigorous standards for biomarker qualification by requiring consistent directional effects across independent cohorts.

Markers testing

Determined R/NR biomarkers were used to calculate log ratios, followed by statistical assessment using the lmer function from the lmerTest package [Kuznetsova et al., 2017] in the R statistical environment. Post-hoc testing was performed using Tukey’s method with the multcomp package [Bretz et al., 2011], applying Bonferroni correction for multiple comparisons. Model residuals were tested for normality using the Shapiro-Wilk test (shapiro.test function in R). Additional testing of marker OGU included investigating the ability of calculated log ratios to separate experimental groups using logistic regression using the method of leave-one group out cross validation. Gene-set enrichment analysis (GSEA) implemented in GSEA function from clusterProfiler package [Xu et al., 2024] using to investigate statistical relationships of OGU marker sets across taxonomy and other attributes included food, body site or pathogens.GSEA analysis results were visualized using GseaVis package [Zhang et al., 2025].

Marker clustering with non-gut databases

The presence of non-gut microbes in the stool microbiome, such as those frequently detected on the oral means, has been associated with a degradation of the gut microbiota [Liao et al., 2024]. This phenomenon frequently results in adverse outcomes, including various illnesses. The identification of non-gut microbes was conducted using published microbial MAGs from body sites (oral, skin, vagina) [Pasolli et al., 2019] and food [Carlino et al., 2024]. The body site MAGs were retrieved from here. A total of 9,412 genomes were selected based on specific criteria, including those from the “Body Site” category, such as the airways, nasal cavity, oral cavity, skin, and vagina. Genomes classified as “Age Category” as either “newborn” or “child” were excluded from the analysis. The selected genomes were then clustered with our reconstructed catalog via dRep [Olm et al., 2017] with 95% nucleotide identity. The genome that exhibited a high degree of similarity with body MAG marks from the same source. Consequently, 293 MAGs were designated as body microbes. The food MAGs were retrieved from here, and the archive cFMD_mags.tar.gz contains the folder cFMD_mags_prok, which contains prokaryotic genomes. These 10,112 metagenome-assembled genomes (MAGs) were subsequently clustered via dRep with 95% nucleotide identity, yielding a total of 983 genomes. The subsequent source identification process was consistent with that employed for body MAGs. Of these, 91 were classified as food microbes.

Functional profiling of marker OGU

Functional profiles of OGUs were generated using MetaCerberus [Figueroa et al., 2024] with annotations from the Carbohydrate-Active EnZymes (CAZy) database [Drula et al., 2022] and the Kyoto Encyclopedia of Genes and Genomes (KEGG) [Kanehisa et al., 2025]. Statistical associations between functional features and experimental groups were assessed using logistic regression. Gene set enrichment analysis (GSEA) was applied to evaluate links between gene sets and immunotherapy response, with results visualized using the ggseavis package.

Additional analysis and visualization

For visualization and statistical analysis, we utilized the following R packages: ggplot2 [Wickham et al., 2016] for graphics creation, ggpubr [Kassambara, 2023] for publication-ready plots, hrbrthemes [Rudis, 2020] for enhanced themes, and RColorBrewer [Neuwirth, 2022] for color palette management. String data manipulation was performed using the stringr package [Wickham, 2022]. The data analysis report was generated with Quarto [Allaire et al., 2022], incorporating dynamic tables via DT [Xie et al., 2023] and enhanced download functionality using downloadthis [Sidi, 2023].

Data availability

In this study, we used open access data from the NCBI/EBI Sequence Read Archives, identified by the following BioProjects accession numbers: PRJNA397906, PRJEB22893, PRJNA399742, PRJNA770295, PRJEB43119, PRJNA762360, PRJNA1011235, PRJNA615114, PRJNA866654, PRJNA494824, PRJEB49516. MAGs catalog assembly pipeline were description in https://github.com/JeniaOle13/Cancer_MAGs. All MAGs sequences were deposited in NCBI under accession PRJNA1196825. Songbird Qurro files for QIIME2 viewer were available at Zenodo. Source code and Quarto report available at https://github.com/JeniaOle13/cancer-biomarkers.

References

Results

Data overview

The analysis incorporated 627 gut metagenomic profiling samples obtained from 11 independent external datasets. Patients were stratified by immunotherapy response into two groups: responders (R group, n = 365; 58.2%) and non-responders (NR group, n = 262; 41.8%). Response assessment followed RECIST 1.1 criteria, with the R group including patients showing complete response (CR), partial response (PR), or stable disease (SD) at 6-month follow-up, while the NR group comprised exclusively progressive disease (PD) cases. The study cohort received various immunotherapy regimens, including anti-PD1, anti-CTLA4, or combination therapies. Cancer type distribution revealed melanoma predominance (n = 456; 72.7%), followed by gastrointestinal (GI) cancers (n = 82; 13.1%), non-small cell lung cancer (n = 15; 2.4%), breast cancer (n = 4; 0.6%), ovarian cancer (n = 2; 0.3%), and other malignancies (n = 68; 10.8%). All samples were collected prior to treatment initiation to evaluate baseline microbiota status.

Figure 1. Global distribution of collected samples (n=627) across countries. Circle size represents sample count per region, while color indicates the proportion of patients responsive (or non-responsive) to immunotherapy.

Quality control

Following quality filtering with fastp, aggregated quality metrics were compiled using MultiQC to summarize the FastQC reports (MultiQC are attached to the report below). Before filtering, approximately ~15bn reads (mean ± SD: 24 ± 15 mln reads per sample) were retained for downstream analysis, with an average of 34 ± 18% of reads per sample removed during quality control (see Table 2 for detailed filtering statistics).

MultiQC report

Taxonomic composition

A total of 3649 OGUs belonging to 12 bacterial phyla were detected in 627 samples. The most numerous of these were Firmicutes (n = 2443), Actinobacteriota (n = 590), Bacteroidota (n = 394), Proteobacteria (n = 140), Desulfobacterota (n = 27) and others (n = 55). OGU abundance profiles of stool samples presented in Table 3.

Marker OGUs discovery

Analysis of marker OGUs using Songbird and Qurro identified 793 significantly differentiated features between R and NR experimental groups. Of these, 326 were enriched in R, while 467 were enriched in NR. The most prominent microbial species enriched in responders included Blautia wexlerae (n = 49), Faecalibacterium prausnitzii (n = 27), Fusicatenibacter saccharivorans (n = 15), Gemmiger qucibialis (n = 14), and Blautia faecis (n = 8). Conversely, NR exhibited enrichment of distinct microbial taxa, including F. prausnitzii (n = 17), Veillonella parvula (n = 7), and Dialister invisus (n = 7). Although F. prausnitzii was detected in both groups, its significantly higher abundance in responders (27 vs. 17 OGUs) suggests potential strain-level or functional variations contributing to treatment outcome. The robust association of Blautia species, particularly B. wexlerae, with positive response highlights its potential as a predictive biomarker. Notably, the enrichment of oral bacteria V. parvula and D. invisus in NR suggests a potential link between oral microbiome composition and immunotherapy efficacy. Detailed descriptions of marker OGUs are provided in Table 4. Songbird coefficients for each OGU marker are presented in Figure 2.

Statistical analysis utilizing linear mixed-effects models, accounting for dataset variables, revealed significant differences between R and NR groups (p < 0.0001). The inclusion of the dataset variable as a fixed effect was statistically significant (p < 0.0001), justifying the employed modeling approach. Residuals were normally distributed (Shapiro-Wilk test, p < 0.0001), and the observed effect size was large (Cohen’s d = 0.91). Log ratio values from the metagenomic samples are detailed in Table 5, with their distribution visualized across experimental groups in Figure 3. A predictive classifier, based on log ratios and logistic regression with Synthetic Minority Oversampling Technique (SMOTE) balancing, demonstrated a receiver operating characteristic (ROC) area under the curve (AUC) of 0.72 ± 0.12 (Figure 4), with other performance metrics detailed in Table 6 and Figure 5.

Saving 8 x 3 in image

Figure 2. The barplot depicts Songbird differential abundance coefficients (MEAN ± SD) for marker Operational Genomic Units (OGUs), with the x-axis representing individual OGUs and the y-axis showing the magnitude and direction of association coefficients.

lmer test results


     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lmer(formula = log_ratio ~ response + (1 | dataset), data = df.log_ratio)

Linear Hypotheses:
            Estimate Std. Error z value Pr(>|z|)    
R - NR == 0   1.7005     0.1508   11.27   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- bonferroni method)

Dataset variable inportance

ANOVA-like table for random-effects: Single term deletions

Model:
log_ratio ~ response + (1 | dataset)
              npar  logLik    AIC    LRT Df Pr(>Chisq)    
<none>           4 -1274.6 2557.1                         
(1 | dataset)    3 -1283.1 2572.3 17.119  1  3.511e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residuals of the model normality testing


    Shapiro-Wilk normality test

data:  residuals(lmer.log_ratio)
W = 0.98021, p-value = 1.817e-07

Effect size estimation

Cohen's d |         95% CI
--------------------------
-0.91     | [-1.08, -0.74]

- Estimated using pooled SD.

Figure 5. Log regression prediction metrics plot.

Marker OGUs features

Gene set enrichment analysis revealed enrichment patterns based on taxonomic features (phylum, genus) and microbial origins – specifically, body sites excluding the intestine (mouth, skin, vagina) and the environment (food). OGU identified with 98% nucleotide identity were grouped by taxonomy or other characteristics to facilitate the study of enrichment within specific groups. According to obtained results, Bacteroidetes phylum was upregulated in R group as well as Proteobacteria in NR group (GSEA adj.p < 0.01). Стоит отметить, что list of core enrichment OGU of Proteobacteria included potential opportunistic species such as Klebsiella quasipneumoniae, Klebsiella michiganensis, Enterobacter ludwigii, Enterobacter kobei, Citrobacter youngae, Citrobacter portucalensis, Citrobacter freundii, Haemophilus sp001815355. At the genus level, Blautia, Bacteroides, Fusicatenibacter, Gemmiger, and Faecalibacterium were enriched in R group as well as Veillonella, Prevotella, Dysosmobacter, Acetatifactor were enriched in NR group (GSEA adj.p < 0.01).

Clustering 3,816 OGUs from this study with 9,412 MAGs derived from different body sites and 10,112 food-derived MAGs revealed that 91 OGUs exhibited complementarity to food-derived microbes, and 293 were complementary to those from non-intestinal sources, at 95% nucleotide similarity, which corresponds to the species level. Based on the data, NR-linked OGUs were enriched in oral and food-derived microbial communities. Core food-derived OGUs comprised Bifidobacterium animalis, C. freundii, E. kobei, Enterococcus faecalis, Enterococcus faecium, Escherichia coli, K. michiganensis, K. quasipneumoniae, Lactobacillus gasseri, Limosilactobacillus reuteri, and Megasphaera elsdenii, showing substantial overlap with core enrichment patterns observed in Proteobacteria. Core enrichment patterns within the oral cavity included OGU representatives of Veillonellales and Enterobacterales, specifically Anaeroglobus micronuciformis, Escherichia coli, Veillonella parvula, and Haemophilus sp001815355.

Functional analysis

GSEA: CAZy

Figure 10. Links between CAZy core enrichment GH and marker OGU.

GSEA: KEGG